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O ■ Abstract 

| The dynamics of a one-dimensional self-gravitating medium, with initial density 

almost uniform is studied. Numerical experiments are performed with ordered and 
with Gaussian random initial conditions. The phase space portraits are shown to be 
qualitatively similar to shock waves, in particular with initial conditions of Brownian 
^ ■ type. The PDF of the mass distribution is investigated. 
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1 Introduction 



In this paper we study the dynamics of a self-gravitating medium, the initial density of 
which is almost homogeneous, and of which the initial velocities of all fluid particles are 
small. 

The study is performed in one spatial dimension where gravitational force is a La- 
grangian invariant. The equations of motion for the density and velocity perturbations 
can therefore be written such that the net acceleration of fluid particles are Lagrangian 
quasi- invariant, if the perturbations are represented as discrete mass points moving on a 
stationary and uniform background. This system of particles is then brought forward in 
discrete time steps from one collision to the next, which, as shown already by Eldridge and 
Feix , permits the construction of an exact (up to round-off errors) and fast numerical 
integration scheme. 

One reason to study the one- dimensional case is, that we wish to study the dynamics 
starting from initial conditions that are Gaussian random fields with power-law spectrum, 
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and the statistical properties of these solutions. By analogy with recent work on the 
adhesion model BO, 133], £TT|, we then expect that there is a wide range of spatial scales in 
the solutions at late times. To get good statistics, while still resolving fully all details of 
the motion, a one-dimensional model is more convenient than simulations in two or three 
dimensions. 

Our motivation to study gravitational dynamics with such random initial conditions 



is that this is as a model of structure- formation in the early universe, see [34, 25, 36]. We 
will return to a discussion of background material in cosmology in section A second 
reason for a detailed study of one-dimensional case is then that, in general position, 
gravitational collapse accentuates asymmetries in the velocity and density fields. The 
first structures to form are blinis or pancakes, thin in one direction and of large extent in 
the two others. What we study can then be pictured as the one- dimensional dynamics of 
very large blinis, oriented parallel to one another and colliding and merging when moving 
under their mutual gravitational attraction. For very long times the tree-dimensional 
nature of the motion is no doubt important, but for some time after the first blinis form, 
the one-dimensional approximation should be appropriate. This point has recently been 



extensively discussed in |[33|| , in the context of the adhesion approximation. The one- 
dimensional model can hence perhaps also give a quantitatively correct description of the 
clustering of mass as function of time as long as we consider the largest structures at 
every moment of time. 

The main new results of this paper are as follows: we rederive the Rouet et al. math- 
ematical model in a way which appear transparent to us, with particular emphasis on 
initially localized perturbations. We discuss the differences between structures observed 
on a uniform and homogeneous background, and those in a finite medium without a back- 
ground. As has been observed previously J^], |2"Bfl , phase space portraits, starting from 
ordered initial conditions, or random initial conditions, consist of smooth one-stream in- 
tervals (with ordered initial conditions), and short intervals with multi-stream solutions 
and high mass concentration. In particular, with Gaussian initial conditions of strong 
spectral support at low wave numbers, e.g. of Brownian type, we find qualitative similar- 
ities to the mass distributions in the adhesion model, i.e. ramps and mass concentrations 
of all sizes at all scales. We try to quantify the mass distribution by computing scaling 
exponents and mass histograms in octaves. 

Recent mathematical investigations of finite self-gravitational systems, without a back- 
ground, are (ordered initial conditions) and (random initial conditions, but of 
another class than we use), for a discussion of applications of this system to structure for- 



mation in the universe, see |12[. For Vlasov-Poisson equation in one dimension (mutual 
repulsive interaction), see f20 |. 

The paper is organized as follows. In section |2| we summarize standard material in 
cosmology, and discuss references to recent observational data on the anisotropies of the 
cosmic microwave background radiation. In section we derive the Rouet et al. solution 
of ID self-gravitating systems, and discuss initial conditions appropriate for our system 
and admissible boundary conditions. In section f| we study the dynamics starting from 
ordered and random initial conditions, and in section |5| we investigate properties of the 
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mass distribution. In section |6| we sum up and discuss our results. 



2 The cosmological and observational setting 

In our immediate neighborhood today, the universe is neither homogeneous nor isotropic. 
Sources of electro-magnetic radiation in any frequency band are distributed in a markedly 
random and clustered manner over the sky. On the other hand, at large scales the universe 
is generally taken to be homogeneous and isotropic. The hypothesis of an early almost 
homogeneous and isotropic state of the universe rests on that it agrees with the whole 
body of theory of the hot Big Bang, and with the observed 3K black-body background 
radiation. It is therefore natural and standard to assume that the structures observed 
today are due to instabilities in an initially almost homogeneous self-gravitating medium 

The study of such instabilities has a long pre-history, going back already to New- 
ton The first quantitative investigation of the instability in a static medium with 



non-zero pressure was performed by Jeans [14], who derived his famous formula that 



perturbations of wave-length larger than Aj = v s ^n/Gp are unstable, where v s is the 
sound speed, G the gravitational constant and p the density. Such perturbations hence 
grow exponentially in time (in the Jeans theory), while perturbations of smaller wave- 
length oscillate and do not grow. The Jeans length can be related to a Jeans mass of 
Mj = ^pAj. In an almost homogeneous gas, regions of increased density with mass 
larger than the Jeans mass collapse gravitationally, while concentrations of smaller mass 
oscillate acoustically. 

While Jeans' analysis immediately applies to gravitational collapse of a finite object, 
it is not complete when referring to a medium of infinite extent. Since a sufficiently large 
mass will always be unstable, the unperturbed state assumed in the Jeans formula, an 
infinite self-gravitating medium with gravitational self-interaction and constant density, 
cannot exist in classical physics. On the other hand, in general relativity the Friedmann 
solutions to the Einstein equations describe unbounded universes with spatially constant 
density. These solutions are given in terms of a cosmological length scale a, which changes 
with cosmological time t. On sufficiently small length scales (much less than a), and on 
sufficiently short time-scales (faster than ") general relativity is well approximated by 
Newtonian gravity, and Jeans' analysis of an infinite medium is hence well-founded in 
this way. 

The linear theory of small perturbations around the Friedmann solutions in general 



relativity was developed by Lifshitz, and is described in detail in Weinberg |34| and 
Zeldovich & Novikov |36| . If we assume a hydrodynamic description of the matter fields, 
the perturbations can be classified as scalar, vector and tensor. The last ones correspond 
to gravitational waves, which have no counterpart in the classical theory, and will be left 
aside in the following. 

Although the relevant calculations are involved, to a considerable degree the results 
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can simply be described by introducing the co-moving wavelength 



?(*) = q(h 



a(i ) 



(1) 



Then the perturbations in density and proper velocity (scalar and vector perturbations) 
grow as 



fMt),t)^e^ x ^ t,)dt 'fMto),t ) (2) 

where Xk{q, t) is the instantaneous growth rate in the Jeans theory of wave- vector q at time 
t, and k labels the linear modes. Equation (Q) does not agree completely with Lifshitz' full 
solution, but is quite close in standard cosmological models. For an extended discussion, 
see @. 

The linear modes can be classified into decaying modes, which are fields of incom- 
pressible proper velocities, and growing modes, which are linear combinations of density 
modes and potential proper velocities. If a mode actually grows or decays at time t also 
depends on whether the co-moving wave-length is smaller or greater than the instanta- 
neous Jeans length. On sufficiently large scales, the perturbations surviving the linear 
regime are coupled density and potential proper velocity fluctuations. 

The outcome of these considerations is that at length scales much smaller then the 
radius of the Universe, but much larger than the Jeans' length, structure formation, as 
long as the solutions are single-stream, is governed by the following system of equations: 



d t p + 3-p + ■ (fi- 



ll 







d t it + -it + -(it ■ ^)lt = — -^ip 



V 2 tp = 4nGa 2 (p - p b ) 



(3) 



If we also assume that initial rotational proper velocity fluctuations have been damped 
out in the linear decay, then it is a potential field at some time to- ^ is the gravitational 
potential generated by the source p — p&. We note that p — Pb can be both negative and 
positive, and is zero on average. On short time scales we can take a constant, and by a 
change of scale we can set it equal to one. On the kinetic level the system is then described 
by the following Jeans- Vlasov-Poisson equation, valid also after caustic formation: 
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(4) 



4ttG(p - p b ) 



The initial conditions of equations (p|) and (H) can be taken to be the fluctuations at 
some stage of linear growth. It is a remarkable fact that the fluctuations at one particular 
time during linear growth are in fact observable in the fluctuations of the 3K blackbody 
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background radiation |23|, |24|, |32|, |TJ|]. At red-shift z ~ 1000 (age of the Universe 10 5 
years), photons fell out of thermal equilibrium with electrons and nuclei; what is observed 
today as 3K black body radiation is the red-shifted spectrum of photons that were in 
equilibrium with matter at that time. COBE observations measure the temperature of 
the blackbody radiation with a beam width of 7°, and detect mean square variations of 
about 30fiK (one part in 10 5 ). Laid out in the sky, COBE can thus be said to distinguish 
a spherical grid of 50 x 50 patches, i.e. a decade and a half in each direction. Experiments 
in the near future (MAP, Planck Surveyor) are expected to increase the angular resolution 
of COBE by more than one order of magnitude. Over the range of COBE, observations 
are in agreement with the Harrison-Zeldovich[|T3|, |35[| prediction of Gaussian initial density 



fluctuations with spectrum |32], [17 



P(k) = Ak n 



n 



1 



(5) 



At scales smaller than about 1°, theoretical arguments predict deviations from fl5|). In cold 
dark matter-dominated models these are determined by the CDM transfer function, which 
at intermediate scales, in the present universe in the range 10 — 50/i~ 1 Mpc, gives a plateau 
where fluctuations are also Gaussian as in (|5|), but with n ~ — 1 ||. We remark that in 
practically all cosmological models, the spectrum (|5|) is not expected to be valid in an 
arbitrary wide range, but to be modified at smaller scales. For further recent discussions 
on the assumed limits of validity of and prospects of experimental observations of such 
deviations, see e.g. |32|, [15], [19|. 



3 The generalized Zeldovich solution: a Lagrangian 
integrable model 

We now turn to a Lagrangian integration scheme where initial mass density p is con- 
centrated on a discrete set of particles. The algorithm which we are going to describe 
was invented in plasma physics (i.e., equation (f|) with the opposite sign of G) in the 
early '60ies, and already previously used in simulations of one-dimensional self-gravitating 
systems |5], |27|, |28|. The derivation we will give stresses localized initial perturbations. 
Boundary conditions are hence those of an unperturbed quiescent state to the left and 
the right of the perturbation. As the perturbation develops it will typically spread, and 
move into the quiescent state, thus inducing that to move. The main ingredient in the 
derivation is a regularisation of the differences between the (formally divergent) forces 
pulling the particle to the left and to the right. Since the physical meaning of an infinite 
classical self-gravitating system is that of an approximation (on sufficient small scales) to 
a system governed by general relativity, and since the speed of propagation of the grav- 
itational interaction in general relativity is finite, localized perturbations are relevant in 
this problem. 

We now take space one-dimensional. The mass density distribution corresponding to 



5 



a system of point-like particles is specified by: 



p(x,t) = ^m k 5{x - x k {t)) 



(6) 



We require that the average density of the point-mass system is the same as the back- 
ground density, i.e. 



Lim/ 



1 

2L 



^2 mk 

k:x k (t)<=[-L,L] 



PO = Pb 



(7) 



and that the perturbation is localized in a weak sense, such that as L tends to infinity 
the measure of the point masses in [L, L + 1] tends to the uniform measure in the same 
interval, with convergence faster than j^. The one-dimensional gravitational potential is 
formally given by 



ip(x, t) = 2G dy\x - y\ (p(y, t) - p b ) {formally) 



(8) 



With the initial conditions under consideration the integral is convergent at infinity. The 
regularisation referred to in the beginning of this section is therefore to take 



ip{x,t) = 2GLim L _ 



mi\x — xi\ 

l:xi£[-L,L] 



Pb 



[x - L) 2 + {x + Lf 



(9) 



where the limit exists by assumption. As now ip is a known and well-defined function of 
position, the equations of motion of the point masses are 



dip{x) 
dx 



-j-k 



-2G 



Lim L _ 



E 

l:xi&[-L,L] 



resign - x{) — 2p b x k 



(10) 



with initial conditions Xk\t=o = x k {0) and Xk\t=o = u {xk{0)). Equation ([10|) expresses 
that the force acting on particle k is equal to the difference of net mass, in excess of the 
background, to the right and to the left of that particle. Hence, were L infinite, it would 
formally be the sum of four terms (two with positive and two with negative signs), each 
of which would be infinite. Our regularisation acts on the net force to the left and to the 
right separately, by the requirement that the initial perturbation is localized. 

With periodic boundary conditions our regularisation would not give a unique answer 
since the limit (|9|) would not exist. One could instead regularise separately the net force 
from particles and net force from the background, for instance by taking the net force of 
each kind initially to be zero, or assuming that the gravitational force actually is of finite 
range, and then take the limit when that range increases. Using ( |I0D necessarily assumes 
implicitly a regularisation, of which ours is a physically transparent one, at the price that 
we consider only sufficiently localized perturbations. 
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We now want to transform fllOD to the Eldridge-Feix scheme. The acceleration acting 
on point mass k at the initial time is 



a(sfc(0),0) 



-2G 



Lim L ^oo 2J m;sign(x fc (0) - Xi(0)) - 2p b x k (0) 

l:xi£[-L,L] 



'ID 



From ([TOj) follows that as long as no other particle overtakes particle k we have the simple 
evolution law 



x k {t) = a(xfc(0), 0) + 4Gpb(xk(t) — x k (0)) (time t before collision) 



(12) 



When particle I overtakes particle k, the gravitational force acting between them changes 
sign. We may therefore write the equations of motion at an arbitrary time t as 



x k (t) = a(x k (0), 0) + AG Pb (x k (t) - x k (0)) + 
-2G^2mi [sign(x fc (t) - xi(t)) - sign(x fc (0) - xi(0))] 



(13) 



If units are chosen such that AGpb = 1 (see below), the mapping from one collision 
to the next is then 

+ ^M±MM exp ^ ( . M)+ (14) 

where x k (t l coll ) and ak(t l coll ) are the velocity and the acceleration felt by the fc'th particle 
at the preceding collision. 



4 Numerical experiments 

A characteristic scale in time is set by the Jeans' frequency, related to the gravitational 
constant G and the mean density p by 

Uj = y^Gio (15) 

In the following we will measure time in units of ujJ 1 . Since the mean density is equal to 
the background density pb, this means that we choose a scale in time such that the value 
of the product AGpb in flllf ) is one. 
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The model of self-gravitating particles without pressure is valid on spatial scales much 
greater than the Jeans length Aj. The intrinsic spatial scale Aj is thus for us zero. Char- 
acteristic spatial scales can hence only come from the initial conditions. If the initial 
perturbation is ordered (smooth function of spatial coordinate), and in analogy with ([!]) 
and (|T0|) has support on an interval [— L, L] then all scales can be measured in terms of L. 
Temporal and spatial scales tuj 1 and L imply however a velocity scale Luj, and we can 
therefore separate ordered initial velocity perturbations as to if they are large or small 
compared to Luj. An initial density perturbation Sp is naturally measured in units of 
the mean density. We choose to normalize mass in terms of the mass M initially involved 
in the perturbation, which means that Pq, and therefore pi,, is The constant G in ( |i~3l) 
is ujj-rjj, in our units hence |. 

The explicit formula (0) allows immediately for a fast numerical scheme with opera- 
tions count O(N) per collision, where iV is the number of particles; this is scheme of 
Eldridge and Feix|7]]. Elsewhere we will discuss an asymptotically more efficient version 
of the algorithm with an operation count of 0(lnN) per collision 0. We note that the 
computer time needed to advance the system up to time t is proportional to the product 
of the number of operations per collision and N co u(t), the number of collisions up to time 
t. If with initially smooth velocity and density perturbations, represented as N discrete 
particles, the mean time between collisions depend on JV as ^, then in the algorithm 
used here the operations count in advancing the system up to intrinsic time t would be 
N 2 . However, due to stretching and separation, the discretization becomes gradually a 
less accurate description of the smooth field to be approximated; this effect is far from 
uniform in phase space, and the putative 0(N 2 ) operation count therefore only holds 
approximately at an initial stage. 

In fig [l] we show an initially smooth velocity perturbation (insert) of sinusoidal shape 
in the interval [—1,1]. The initial density perturbation is zero. The full system thus 
consists of a uniform stationary background with density p&, and a number of particles 
as in (HD, distributed uniformly with respect to the same density, and at rest outside the 
interval. As long as no particle from the inside has reached the boundaries all three terms 
on the right-hand side of equation ([13]) remain unchanged for a particle outside [—1, 1], 
and these particles thus remain at rest. The main figure in fig [I] show the solution at time 
4.0 uj 1 , at a resolution where obviously the continuum description is still applicable, and 
where no particle from inside has reached the boundaries. In phase space the distribution 
has support on a curve of spiral shape. If initial velocity is large or compared to Lujj, the 
fastest particles typically reach the boundaries before a spiral is formed, while if initial 
velocity is much smaller than Lujj one or more turns of the spiral form before any particles 
reach the boundaries. These observations have already been made by Rouet et al , 
and show a qualitative difference in the solutions depending on the initial kinetic energy. 

In fig ^ we show the development of the system in fig [I] at successive later stages. The 
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most distinct features are the formation of a high-density cluster in the middle, and two 
expanding high-density fronts to the left and right, but with very low density in between. 
The two fronts form when particles from the interval collide with particles that were ini- 
tially at rest. The first particles to do so overtake the same mass in resting particles and 
background, and will therefore feel no change in the acceleration. The resting particles 
that have been overtaken will on the other hand feel an acceleration towards the front 
that has passed. They will hence start moving towards the front, making non- leading 
particles in the front overtake more background than particles, and therefore also feel 
an increased acceleration. Both these effects lead to a high density concentration at the 
front, large force gradients, and strong stretching in phase space. In fact, it is clear that 
with the resolution used in fig. ||, at the fronts the continuum description is already lost. 

In the center of the middle cluster the density is several times the background, and 
the spirals are round, as in dynamics without the background, while in the outer parts, 
where 5p is comparable to pb, the spirals are deformed, as above on fig [l]. 



We now turn to random conditions. We made three different choices of velocity as function 
of spatial position: Brownian motion; fractional Brownian motion with Hurst exponent 
equal to zero, and white noise. All three are random Gaussian fields with power law spec- 
tra as in eq. of which we choose white noise and Brownian motion because they are 
Markov processes, and have been investigated in the context of the adhesion approxima- 
tion [30, plj| , white noise also as a reference case, because it has earlier been investigated 



by Rouet and collaborators |27, 28) , and the fractional Brownian motion as an interesting 
intermediate case. 



The translation between exponents in (|5]) (density perturbations in 3D) and our initial 
conditions (velocity perturbations in ID) is as follows: a Gaussian random function / with 
stationary increments has second order correlation function < \f(x + I) — f(x)\ 2 >~ l 2h , 
where the Hurst exponent h £ [0, 1[ is related to the scaling exponent n of the spectrum 
E(k) ~ k n , through 2h + n + (D — 1) = —1. A scaling exponent n\D in ID is therefore, for 
these purposes, analogous to a scaling exponent n 3 o = uid — 2 in 3D. A density perturba- 
tion Spk is in the linear regime on the other hand tied to a velocity perturbation Vk ~ kpk- 
We therefore have = + 2, and combining the two relations nf^ = n^. Our 
intermediate case with Hurst exponent equal to zero hence corresponds to the Bardeen- 
Bond-Kaiser-Szalay intermediate spectrum with nf^ = — 1, and will in the following be 
referred to as BBKS initial conditions. White noise and Brownian motion correspond to 

nil = and n 3D = -2. 

Gaussian random fields are conveniently generated in the Fourier space representation. 
At time t = particles are uniformly distributed on an interval of size I with unit spacing 
Ax — i/N and velocity 

v(x) = ^ k v k e ihx (16) 

The sum over k extends from — ir / Ax to n / Ax in steps of 2 7r / i and v = v l- The 
Fourier components of positive k are then chosen as a random Gaussian independent 
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Figure 1: Velocity field vs. position at time t = 4.0 ujj . Number of particles is iV = 
15 x 10 3 , initial sinusoidal conditions in lower right corner insert. 
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Figure 2: Velocity field vs. position of the system in fig. |T| at a later time, well into the 
non- linear regime (t = 11.08 uj 1 ). The insert is a blow-up of the central region. Three 
high-density peaks are clearly developed. 
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variables with variances: 

U-2h-l 

< \ Vk \ 2 >= — (17) 

1 1 2 Ax K J 

The field generated by (|T6|) and (|17|) will be periodic with period i. If /i is in the interval 
[0, 1[ the field will have stationary increments on length scales much less than £, and on 
these scales approximate a Fractional Brownian motion with Hurst exponent h, while if h 
is in the interval ] — 1, 0[ the field itself will be stationary on length scales much less than 
£, and approximate the derivative of a Fractional Brownian motion. We choose i = 2L, 
so that particles are initially distributed in a box [-L, L}. 

We will now discuss units of time and space with random initial conditions. One scale 
is 2L, another is the ultra-violet cut-off luvi which in our case is at least as large as the 
initial particle distance Ax, and a third is an infra-red cut-off lm, which is not larger than 
2L. With h in the interval [0, 1[ we have in law 

v(x + /) - v(x) ~ v uv (-^-) h (18) 

l-uv 

where vuv is the size of typical velocity fluctuations on the ultra-violet cut-off scale. The 
typical overturning time at scale / is then 

U ~ l -™(-Ly~ h he[o,i[ (19) 

vuv L uv 

The time can be measured in units of inverse Jeans' frequency, and be small or large 



in those units. Equation (|19f) then predicts that the characteristic time to form a struc- 
ture of size I increases with I, such that small scales form first. At sufficiently large scales 
I the initial velocity fluctuations will be small compared to luj, and we hence expect to 
see the central spiral structure of figs |I] and but not much of the fronts. 

On the other hand, if the spectral exponent is larger than —1 (Hurst exponent less than 
zero), then the initial velocity field is homogeneous, characterized by a rms velocity vrms, 
which is on the order of vuvi an d we therefore expect the simpler result t\ ~ (l^) - 
At scales much larger than lyy velocity fluctuations are again much smaller than luj, 
and we therefore expect mainly to see the central spiral structures of fig [l] and fig |2|. We 
remark that if we reason by analogy, and assimilate these structures to shock waves in the 
adhesion model, which trap fast particles, then the characteristic times will be longer, and 
in fact again of the form (p!9|). The turn-over time ti then grows faster then linearly with 



I [fll, [1C[| . Unfortunately our present resolution is not sufficient for a precise determination 
the typical size as function of time and of the temporal development of the spectral shape 
of the perturbations. 

In figs. |3|, [|, |5|, ^ we show the phase space portraits with Brownian motion, white noise 
and BBKS initial conditions, respectively. As expected, we see many spiral structures, 
small and large, and fronts at the left and right boundaries, except in fig. f|, included as a 



11 



consistency check (see below). In qualitative agreement with the adhesion model, we also 
see "ramps" with low density, and where velocity is an increasing function of position, 
interspersed with regions of high density. By visual inspection of an agglomeration of a 
certain scale, it appears that the velocity to the right of such an agglomeration is less 
than on the left, such that velocity has negative gradient through the agglomeration, just 
as in shocks in the adhesion model. This picture is however complicated by the fact that 
there are agglomerations of different sizes, and that they are not simply ordered from 
left to right. A certain stretch of phase space, say just to the left of some structure, is 
on some other scale included inside a much larger structure; this effect appears to be 
especially pronounced with Brownian motion initial conditions. In the adhesion model 
such microstructures are of course hidden inside the shocks. 



The difference between figs. |] and [5] is that in fig. [5] we include, following our gen- 
eral approach, a medium of quiescient particles to the right and left. As in fig. ^ we 
then find two expanding fronts, where the particles of initial velocities around Vrms i n 
the interval collide with the particles initially at rest. To show that the fronts do not 
qualitatively change the dynamics in the central region, we show for comparison in fig. |] 
a system with only background outside the interval. Particles that escape the central 
interval are then accelerated outwards by the anti-harmonic force in (0), and give phase 
space portraits where velocity v is an increasing function of position x outside the interval. 



Another way to eliminate the fronts would be to take a perturbation which is similar 
to white noise, but of which the envelope would go smoothly to zero outside an interval of 
length 2L, which can be achieved with essentially flat spectrum at wave numbers larger 
than For Brownian and BBKS initial conditions we do not have to worry about the 
fronts. Without changing any many-point statistics we can constrain a realization us- 
ing flTE| ) and (O) to vanish at one point, say x = —L, and then by periodicity also at 
x = L. By (|T8|) the typical velocities in the center of the interval [— L, L] will then be 
much larger than at the boundaries, and the most prominent structures are found there. 

To quantify how accurate is a finite-particle description we looked at PDF of the 
inter-particle distances at different times. In fig. ^ we show the PDF of the inter-particle 
distances to the power one-half (insert on the right). In each case a peak is clearly 
observed: the positions of these peaks are independent of the initial realizations and 
evolution time. The main plot in fig. |7] shows the PDF of the rescaled quantity T = (|f ) 1 / 2 
where the constant (3 = 0.6 was obtained by a numerical fit. We have no good explanation 
of the observed scaling behaviour at this point. 



5 Mass analysis 

Qualitatively, the velocity field recalls the results obtained in the framework of the ad- 
hesion model |TI| |3(| |33fl , where "ramps" appear that separate shock regions where fully 
inelastic collisions among the Lagrangian particles occur. It is interesting to try to in- 
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Figure 3: Velocity field vs. position at time 4.76 uij , starting from single-speed Brownian 
motion initial conditions (small upper right insert). Number of simulated particles is 
N = 8192. 




Figure 4: Velocity field vs. position at time 7.37 ooj , starting from single-speed white 
noise initial conditions (lower insert). The model is modified, as described in the main text, 
to contain only background and no particles outside the interval of the initial perturbation. 
Number of simulated particles is N = 8192. 
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Figure 5: Velocity field vs. position at time 7.37 ooj , starting from single-speed white 
noise initial conditions as in fig. |] inside the interval of the perturbation, and quiescient 
particles to the right and left thereof. Note the formation of fronts as in fig. |2|. Number 
of particles initially involved in the perturbation, with total mass one, is N = 8192, and 
the total number of simulated particles, including the ones to the left and right, is 16384. 




Figure 6: Velocity field vs. position at time t = Q.7ujj starting from single-speed BBKS 
initial conditions (upper right insert). Number of simulated particles N = 8192. 
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troduce quantitative characteristics of the mass distribution in order to make a closer 
comparison. 

One main result of [30| and ||, is that for h £ [0, 1[ the inverse Lagrangian map, initial 



versus final position, has a bifractal structure, similar to that of the Devil's staircase 
construction using the standard triadic Cantor set. Except for a set of measure zero, all 
Lagrangian initial points land on shocks, which in Eulerian coordinates are at definite 
shock locations, and where therefore all of the mass is concentrated. The number density 
per unit length of shock locations holding mass m is in a well-defined range distributed 
as a power-law m~ h ~ x . Most mass therefore lies in the largest shocks formed at a given 
time, but the number of smaller shocks is divergent. For the Brownian motion initial 
condition it can be proven, and for the other cases it is conjectured, that Eulerian shock 
locations are almost surely dense. Most of these shocks are still small. However, since 
all mass initially uniformly distributed after an arbitrarily short time falls into a shock, 
the mass contained in an interval after a finite time is in this model only made up of the 
mass in the shocks actually inside the interval. The mass measure of the adhesion model 
with these initial conditions is therefore bifractal, which can be quantified by the scaling 
exponents of its moments 

M(q,l) = Z? =1 p(xi) q l~r {q) (20) 

where / is the length of the coarse-graining mesh, B is the number of boxes in the mesh, 
and L = Bl. The sum is normalized such that r(l) = 1 and r(0) = 0. At sufficiently 
large q, where the threshold lies at h, the sum in ( p0|) is dominated by a small number 
of terms, corresponding to the intervals containing the larger shocks. The exponent r 
in this range is then one. At small q (^0]) would instead be dominated by almost empty 



intervals, each of which carries a mass , and of which there would be Z -1 in number. 
The scaling exponent in this range would hence be r(q) = |. With n in the range [—1,1] 
(which formally corresponds to h in the range [—1,0]), of which one case is white noise 
{n — 0, h = — |) all the above statements remain true in the adhesion model, but some- 
what trivial. Shock locations holding mass m are still distrubuted as m _/l_1 , but since h 
is now negative, most shocks are within an octave in size of the largest. Shocks are not 
dense, the mass distribution is almost surely concentrated on a finite number of points 
per unit length, and r(q) = 1 for all positive q. 



If the bifractal scaling behaviour of the moments in (|20| ) would be observed also for a 
self-gravitating system, then the two models would in this sense be equivalent. Possible de- 
viations from bifractality would on the one hand be intrinsic effects of the self- gravitating 
dynamics. In fig. ||| the local scaling exponent r = r(q, I) is shown as a function of I. One 
problem, discussed at length in for the adhesion model, is that with a finite number 
of particles, true scaling behaviour can only be observed in a range where most intervals 
actually contain a particle. At smaller mesh sizes, most boxes will be empty, and ( pOD 
would be dominated by a small fraction at any positive value of q, which gives the spurious 
result r(q) = 1, also for q in the interval [0, h]. The predicted cut-off occurs at l sp ~ e h , 
where e is the numerical mesh size in a simulation of Burgers' equation, and which we 
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in analogy in our case can take to be 2L/N. The value of l sp is thus unfortunately not 
very small. In the present problem mass is not so concentrated, and we could expect to 
perhaps see a slightly wider range. Nevertheless, it is clear that in the range q less than 
one, the only cleanly observed scaling behaviour in I is the spurious one at small I, in our 
case I less than about 10 -4 . 



A more direct quantification of the mass distribution is the mass octave function (MOF), 
i.e. the probability to find a non zero contribution to the mass density as function of 
the mass itself, coarse-grained in octaves. As discussed above, in the adhesion model the 
number of intervals with mass in an octave around m scales as mT h 
behaviour would be observable in all three models we study. 



30, 33 , and this 



At a less detailed level, the MOF measures the degree of dishomogeneity of a proba- 
bility distribution, since in terms of the MOF a uniform distribution corresponds to a 
logarithmic histogram with only one non-zero entry. 



Figs. H [TO] and 11 show the MOF for Brownian motion, BBKS and white noise ini- 
tial conditions. As we see, the observation of a non-uniform density distribution is clearly 
borne out, since the MOF diagrams have support over several octaves in all three cases. 
More interestingly, although statistics and the range are not very large, for white noise 
initial conditions the most massive boxes are the most frequent ones, while for Brownian 
motion less massive boxes are at least equally frequent, in qualitative agreement with the 
adhesion model. The BBKS model seems to be intermediate. 



6 Discussion 

Observed structures in our universe are believed to have been caused by gravitational 
instabilities in an initially almost homogeneous medium. The temperature fluctuations 
of the Cosmic Microwave Background radiation give direct observational access on the 
primordial density and proper velocity fluctuations which seeded these structures, albeit 
as of today only on a limited range of scales. 

The general problem of structure formation by nonlinear deterministic evolution equa- 
tions, acting on random initial conditions, is a topic of much current interest, particularly 
in the context of Burgers' turbulence, where progress has been made using the special inte- 
grability properties of Burgers' equation. In this paper we have investigated the dynamics 
of a one- dimensional self-gravitating medium, as a more accurate model of structure for- 
mation in the universe. Therefore, it is of interest to understand how different are the 
solutions to Burgers' equation and self-gravitating systems starting from the same initial 
data. 

The problem addressed is hence not to determine if these two models are closely similar 



16 



0.08 




0.01 0.02 0.03 

r 



Figure 7: Insert on the right: the PDF of the inter-particle distances to the power 1/2, 
for Brownian initial velocities (squares, t = 6AUJ 1 ), BBKS (circles, t = 7.5UJ 1 ) and white 
noise (triangles up, t = 8.1UJ 1 ). Main figure: PDF of T, where T = (# ) 1/2 and (3 = 0.6. 
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Figure 8: Local scaling exponent obtained as the logarithmic derivative of M(q, I) versus 
/ at time t = 5.8 uj 1 with Brownian initial velocities. Number of particles is iV = 8192. 
The exponent q ranges from 0.125 to 1.875 in steps of 0.125. 
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Figure 9: Normalized mass distribution in octaves for Brownian initial velocities at 

t = 6.0 ujJ 1 . The number of particles is 8192, number of independent realizations 30. 

-l 



m intrinsic 



Each realization was further averaged over 10 4 collisions, about 2 x 10 _3 u; 
time. The height of a column is the fraction of a total number of bins containing mass in 
the range [m, 2m]. The bin size is I = 0.003125, a uniform density would hence correspond 
to a single column at abscissa 0.003125 and ordinate 1.0. 



Figure 10: Normalized mass distribution in octaves for BBKS initial velocities at t = 
5.9 ujJ 1 , Number of particles and bin sizes as in fig. || number of independent realizations 
20. 
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in general, but if they are similar with specific initial conditions suggested by cosmology. 
We have focused on Gaussian random fields with scaling spectrum, of the type originally 
suggested by Harrison and Zeldovich, but with a variable scaling exponent. The main 
obstacle to a detailed comparison is then to solve numerically the self-gravitating system, 
since Burgers' equation is solved directly with the Hopf-Cole transformation. In one spa- 
tial dimension we have been able to use an efficient numerical algorithm which exploits 
a special Lagrangian quasi-invariance of the gravitational force between particle collisions. 

A self-consistent formulation of an infinite self-gravitating system demands a background 
term from the average mass density. This somewhat trivial term changes the properties 
of the solutions qualitatively compared to those of a finite mass concentration, with zero 
background. A finite self-gravitating system has a definite center of mass, and particles 
which are furthest from the center feel the strongest attraction, while in the self-gravitating 
system with background mass far from the initial perturbation feel no attraction at all. 
As soon as the system with background has developed structures of much higher density 
than the background their further development is however quite similar to a finite self- 
gravitating system: this shows up in the formation of a central body of enhanced density 
and spiral phase-space structures. A mathematical difference between perturbations in a 
finite and an infinite self-gravitating system is an effective anti-harmonic term in the lat- 
ter, which appears when particle density thins out. Since if in the infinite self-gravitating 
system density can only be less than the average locally if it is higher elsewhere, a thin- 
ning out assumes an attractive agglomeration at some other nearby position, which can 
be seen as the cause of that extra force. 



One prediction with analogy with the adhesion model is that for initial perturbations 
with strong support at low wave numbers, as our case Brownian motion, we expect to 
see mass agglomerations of very different sizes, while for white noise initial conditions 
we expect to find most mass agglomerations of similar size. Figs. || 
behaviour, while fig. ITU] is intermediate. 



11 indeed show this 



Summing up, we have shown that one-dimensional self-gravitating dynamics can be in- 
vestigated quantitatively in systems with a large scaling range in the initial conditions. 
Nevertheless, the problem remains computationally more demanding than e.g. Burgers' 
equation, and our results on the fractal properties of the mass distribution are not conclu- 



sive. An improved numerical procedure has been proposed by A. Noullez [22], which would 



allow us to significantly enhance the resolution and the statistics. Results of this second 
stage of the investigation will be presented in a forthcoming separate contribution p|. 
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Figure 11: Normalized mass distribution in octaves for white noise initial velocities 
at t — 7.0 cjJ 1 . Number of particles and bin sizes as in fig. || number of independent 
realizations 10. 
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